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1.  Introduction 

One  of  standard  approaches  to  numerical  treatment  of  boundary  value  problems  for  elliptic 
partial  differential  equations  (PDEs)  calls  for  converting  them  into  second  kind  integral  equations 
(SKIEs)  with  subsequent  discretization  of  the  latter  via  appropriate  quadrature  formulae.  Dis¬ 
cretization  of  the  resulting  SKIEs  usually  leads  to  dense  luge-scale  systems  of  linear  algebraic 
equations,  which  are  in  turn  solved  by  means  of  some  iterative  technique,  such  as  Generalized  Con¬ 
jugate  Residual  algorithm  (see  [11,  21]).  Most  iterative  schemes  for  the  solution  of  linear  systems 
of  this  type  require  application  of  the  matrix  of  the  system  to  a  sequence  of  recursively  generated 
vectors.  Applying  a  dense  matrix  to  a  vector  is  an  order  n2  procedure,  where  n  is  the  dimension  of 
the  matrix,  which  in  this  case  is  equal  to  the  number  of  nodes  in  the  discretization  of  the  domain 
of  the  integral  equation.  As  a  result,  the  whole  process  is  at  least  of  the  order  n2,  and  for  many 
large  scale  problems,  this  estimate  is  prohibitively  large. 

In  the  present  paper,  we  describe  an  algorithm  for  rapid  application  of  matrices  resulting  from 
discretization  of  integral  equations  of  scattering  theory  in  two  dimensions  to  arbitrary  vectors.  The 
algorithm  requires  an  amount  of  work  proportional  to  n4/3  where  n  is  the  number  of  nodes  in 
the  discretization  of  the  boundary  of  the  scatterer,  and  when  it  is  combined  with  a  Generalized 
Conjugate  Residual  type  algorithm,  the  resulting  process  takes  very  few  iterations  to  converge, 
leading  to  an  order  n4/3  algorithm  for  the  solution  of  the  original  scattering  problem. 

Reduction  of  boundary  value  problems  for  elliptic  PDEs  to  second  kind  integral  equations  is 
discussed  in  detail  in  [6, 15, 17].  Numerical  treatment  of  SKIEs  in  the  general  case  can  be  found,  for 
example  in  [3],  and  numerical  solution  of  acoustic  scattering  problems  in  two  dimensions  by  means 
of  SKIEs  is  discussed  in  [16].  We  present  an  algorithm  for  rapid  solution  of  integral  equations  of 
classical  potential  theory  (Laplace’s  equation)  in  [18],  and  the  algorithm  of  the  present  paper  can 
be  viewed  as  an  extension  of  the  approach  of  [18]  to  the  case  of  the  Helmholtz  equation.  However, 
the  analytical  apparatus  of  the  present  paper  is  considerably  more  complicated  than  the  analytical 
apparatus  of  [18],  reflecting  the  difference  between  the  behavior  of  solutions  of  the  Helmholtz 
equation  and  that  of  harmonic  functions. 

Remark  1.1 

While  the  algorithm  of  the  present  paper  has  an  asymptotic  CPU  time  estimate  n4/3  ,  it  can 
be  easily  modified  into  an  order  nlog(n)  algorithm  (see  Subsection  4.4).  However,  it  appears  that 
this  modification  would  not  lead  to  significant  improvement  in  actual  calculation  times  for  most 
problems  of  practicable  size  (n  <  20000) . 

2.  Background  Information 
2.1.  Notation 

We  will  be  considering  the  situation  depicted  in  Figure  1.  A  fluid  scatterer  of  arbitrary  shape 
is  imbedded  in  a  two-dimensional  fluid  space.  The  boundary  of  the  scatterer  parametrized  by  its 
length  will  be  denoted  by  7  so  that  7  :  [0,  L)  — ►  i?2  is  a  Jordan  curve,  and  the  image  of  7  will  be 
denoted  by  T  .  The  open  interior  of  7  will  be  denoted  by  Q,  so  that  T  =  dtl.  We  will  assume  that  7 
is  at  least  c2,  i.e.,  that  at  each  point  it  has  at  least  two  continous  derivatives.  The  interior  normal 
to  7  at  the  point  x  =  7 (t)  will  be  denoted  by  N(t).  and  it  will  always  be  assumed  that  [| AT (t)||  =  1. 
The  density  of  the  scatterer  will  be  denoted  by  and  the  speed  of  sound  in  it  will  be  denoted  by 
c,n .  The  density  of  the  containing  space  will  be  denoted  by  pout ,  and  the  speed  of  sound  in  it  will 
be  denoted  by  cout.  We  will  denote  the  angular  frequency  of  the  source  by  ui,  and  its  location  by  xt. 
Finally,  we  will  denote  the  Helmholtz  coefficients  inside  and  outside  the  scatterer  by  and  kout 
respectively  (  as  is  well  known,  k,„  =  w/cm,  koui  =  u >/cout  ).  In  the  notation  introduced  above,  we 
assume  that  pm  and  pout  are  positive  real  numbers,  and  that  cm,  coul  and  w  are  complex  numbers 
such  that  Re(cin)  >  0,  Re(cout)  >  0,  f?e(w)  >  0  and  Jm(w)  >  0,  Im(cin  <  0),  Im(cout)  <  0. 
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2.2.  Single  and  double  layer  potentials. 

For  a  Helmholtz  equation 

+  =  0  (2.1) 

we  will  define  the  field  :  R2  \  {x0}  -*  C1  of  a  unit  charge  located  at  the  point  xo  €  R2  by  the 
formula 

<(*)  =  ffo(*|j*-*0||).  (2-2) 

We  will  define  the  field  4>kXo^  of  a  unity  dipole  located  at  xq  and  oriented  in  the  direction  h  €  R2 
by  the  formula 

<*(»)- -*(t||.-«||).^2jjp.  (2.3) 

For  a  continous  function  a  :  (0,  L]  — *  C1,  the  potential  of  a  single  layer  of  density  a  on  a  curve  7  is 
a  mapping  Pfc, :  R2  — *•  C1  defined  by  the  formula 

rL 

/  4>k{t)(x)o(t)dt  (2.4) 

J  0 

and  the  potential  of  a  double  layer  of  density  a  on  a  curve  7  is  a  mapping  Pj^  :  R2  —*  C1  defined 
by  the  formula 

[L 

Plrix)  =  /  ^{t)N{t)[x)(T(t)dt.  (2.5) 

JO 

Remark  2.1 

Note  that  while  both  It  and  P]^  are  defined  on  all  of  R2,  neither  Pj^  nor  the  derivatives  of 
Pto  are  continous  in  the  neighborhood  of  I\  The  exact  nature  of  their  singularities  is  crucial  for 
the  derivation  of  the  equations  2.13,  2.14  of  the  following  subsection,  and  it  is  discussed  in  great 
detail  in  [16],  [15]. 

2.3.  Acoustic  scattering  in  two  dimensions. 

In  the  present  paper,  we  will  be  considering  the  following  problem: 

For  a  pair  of  continous  functions  /,  g  :  T  — ►  C1,  find  two  mappings  <p  :  Cl —*  Cl ,  \p  :  R2\Cl  — ►  C1 


such  that 

a. 

V2*  +  k2n<t>  =  0  on  0 

(2.6) 

b. 

VV  +  k2ulip  =  0  on  R2  \  Cl 

(2.7) 

c. 

pin  ■<p  +  pout  -rp  =  J  on  r 

(2.8) 

d. 

Q 

-(«!>-  ip)  =  g  on  T 

(2.9) 

e.  ip  satisfies  the  radiation  condition  at  oc.  i.e..  for  any  x  €  R2.  there  exists  c  e  Cl  such  that 


Kin  xp(t  •  x)  •  y/(t)  =  c  (2.10) 

with  k  =  koxu  . 

The  above  five  equations  describe  acoustic  scattering  from  a  two-dimensional  fluid  inclusion  in 
a  fluid  space  in  the  frequency  domain,  and  have  been  studied  in  great  detail  (see,  for  example,  [2], 


[5],  [10],  [16].  Their  numerical  solution,  however,  presents  a  number  of  serious  difficulties,  especially 
for  large*  scale  problems.  Here  we  will  follow  the  approach  of  [16],  which  calls  for  reducing  these 
equations  to  second  kind  integral  equations  and  solving  the  latter  numerically.  As  is  shown  in  [16], 
by  introducing  two  new  unknown  functions  <r,i?  :  [0,1,]  — ►  C1  and  representing  the  functions  <f>,  tp 
in  the  form 

*  =  +  (2.11) 
*-  ■psF&..i+pL,'’  (2.i2) 

the  equations  (2.l)-(2.5)  are  reduced  to  a  pair  of  second  kind  integral  equations  on  the  boundary 
of  the  scatterer: 

-2 »V“  +  0»  +  (/£.„  -PL,)  +  (e°*PL.,  -  e‘"PL°)  =  /.  (2.12) 

»<;=?  +  pi  +  p(pL,'PL.')  +  pipPl-.,  -  p*L,)  -  »■  (2.14) 

2.4.  Iterative  solution  of  second  kind  integral  equations. 

The  system  of  equations  (2.13),  (2.14)  satisfies  the  conditions  of  the  Fredholm  theorems  and 
can  be  efficiently  solved  by  means  of  Generalized  Conjugate  Residual  type  iterative  algorithms  (see 
[11],  [16],  [21]).  Iterative  solution  of  integral  equations  usually  involves  application  of  the  integral 
operator  in  the  left-hand  side  of  the  equation  to  a  sequence  of  recursively  generated  functions. 
Applying  an  integral  operator  to  a  function  numerically  is,  generally  speaking,  an  order  n2  pro¬ 
cedure,  where  n  is  the  number  of  nodes  in  the  discretization  of  the  domain  of  the  operator.  The 
resulting  CPU  time  estimate  for  the  solution  of  the  original  scattering  problem  is  also  of  the  order 
n2  (see  [3]  [16]  ),  which  can  be  prohibitively  expensive  for  large-scale  problems.  The  rest  of  this 
paper  is  devoted  to  constructing  an  algorithm  for  numerically  applying  the  integral  operators  in 
the  left-hand  side  of  the  equations  (2.13),  (2.14)  to  arbitrary  functions  in  a  “fast”  manner,  i.e.  for 
a  cost  less  than  n2  (the  particular  algorithm  we  have  tested  has  an  asymptotic  CPU  time  estimate 
n4/ 3  ). 

Remark  2.2 

Evaluating  integral  operators  in  the  left-hand  sides  of  equations  (2.13),  (2.14)  numerically 
can  be  viewed  as  evaluating  the  fields  and  normal  derivatives  of  the  fields  created  on  the  curve 
•7  by  charge  and  dipole  distributions  on  that  same  curve.  In  the  following  section,  we  develop  an 
analytical  apparatus  for  rapid  evaluation  of  fields  (and  derivatives  of  the  fields)  of  distributions  of 
charges  and  dipoles,  and  in  Section  IV,  we  use  this  apparatus  to  construct  an  actual  algorithm  for 
rapid  evaluation  of  integral  operators  of  the  forms  (2.13),  (2.14). 

2.5.  Asymptotic  behavior  of  Bessel  functions. 

In  agreement  with  the  standard  practice,  we  will  denote  by  Jm  the  Bessel  function  of  the  first 
kind  of  order  m,  and  by  Hm,  the  Hankel  function  of  order  m.  As  is  well  known  (see,  for  example 
[20]),  Jm  are  analytic  on  the  whole  complex  plane  for  all  values  of  m,  while  Hm  have  a  branch  cut 
along  the  negative  real  axis,  and  become  infinite  at  the  origin.  The  asymptotic  behaviour  of  the 
functions  Jm ,  Hm  for  large  m  is  given  by  the  formulae 

lim  Jm{=)  •  (— r  •  sn&n)  =  1  (2.15) 

HI— *00  CZ 


lim  H, 


m— oo 


y/{irm)  _ 

v/2  " 


and 


(2.16) 
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(see  [1],  9.3.1,  9.3.2,  9.1.3).  For  large  z  and  fixed  m,  the  asymptotic  behavior  of  Jm(z),  Sm(z)  is 
given  by  the  formulae 


~  \f{^)  cos{z-^y~^)  =  0{ 


g/TO(*) 

kT 

V~zHm(z)  -  =  O(^) 


). 


when  z  — *•  oo,  as  long  as  Im(z)  >  0  (see  [lj,  9.2.5,  9.2.7). 


(2.17) 

(2.18) 


3.  Rapid  Evaluation  of  Radiation  Fields. 

3.1.  Partial  wave  expansions  of  radiation  fields. 

If  a  function  <j> :  R2  —*  C 1  satisfies  the  Helmholtz  equation  (2.1)  in  an  open  disk  D  of  radius 
R  with  the  center  at  the  point  xo  €  R2  then  there  exists  a  unique  sequence  a  =  (am),  m  = 
0,  ±1,  ±2,  •  ■  ■  such  that  for  any  x  €  D, 


+00 

*(*)-  E  amJm(kp)eimt.  (3.1) 

m«-oo 

In  the  above  formula,  p—\\x  —  io||  and  6  is  the  angle  between  the  vector  x  —  xo  and  the  x  axis. 

If  a  function  rp  satisfies  the  equation  (2.1)  outside  D  and  the  radiation  condition  (2.10)  at  oo 
then  there  exists  a  unique  sequence  f3  —  {/9m},  m  =  0,  ±1,±2,  •  •  •  such  that  for  any  x  €  R2\D, 

+00 

*(*)-  E  0mHm(kp)eiml.  (3.2) 

m«-» 


A  derivation  of  the  formulae  (3.1),  (3.2)  can  be  found,  for  example,  in  [15],  and  we  will  refer  to 
functions  satisfying  the  Helmholtz  equation  as  radiation  fields,  to  expansions  of  forms  (3.1),  (3.2) 
as  J-expansions  and  H-expansions  respectively,  and  to  the  point  xq  as  the  center  of  the  expansions 
(3.1),  (3.2). 

The  following  lemma  is  a  direct  consequence  of  the  formulae  (2.15),  (2.16).  It  establishes  the 
convergence  rates  of  the  expansions  (3.1),  (3.2). 

Lemma  3.1. 

If  Di  C  D  is  a  disk  of  radius  R\  <  R  with  the  center  at  xo  then  there  exists  c  >  0  such  that 
for  any  x  €  Di  and  N  >  |&|  •  Ri, 


W*)-  E  <«(§■)"■ 


(3.3) 


mm—N 


If  D2  D  D  is  a  disk  of  radius  R2  >  R  with  the  center  at  xo  then  there  exists  c  >  0  such  that  for 
any  x  €  R2  \  £>%  and  N  >  |fc|  •  R, 


iv-w  -  £  p„Hm(kPwn  < 

mm—N  1 


(3.4) 


Remark  3.1 
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In  numerical  calculations,  expansions  (3.1),  (3.2)  are  truncated  after  a  finite  number  of  terms, 
and  the  resulting  expressions  are  viewed  as  approximations  to  the  fields  <f> ,  0.  If  we  want  to 
approximate  <f>  by  an  expansion  of  the  form  (3.3)  with  an  accuracy  e  then  according  to  the  above 
lemma,  we  have  to  choose 


N  >  max(|k|  •  Ri, 


-ln(f)  -hln(c) 


)■ 


In (R)  -  ln(f?!r 

Since  logarithm  is  a  very  slowly  growing  function,  for  medium  and  large  scale  problems, 


(3.5) 


max(f?i  -  |fc|, 


-  ln(e)  +  ln(c) 
ln(J2)  -  ln(JZj) 


)  ~  Ri  ■  |fc|. 


i.e.  the  number  of  terms  in  the  approximation  is  almost  independent  of  e,  and  must  be  roughly 
equal  to  |Ar|i2i.  A  similar  calculation  shows  that  for  medium  to  large  scale  problems,  the  expansion 
(3.2)  can  be  truncated  after  approximately  N  =  |fc|i?  terms. 


3.2.  Translation  operators  for  H  and  J-expansions. 

For  a  real  number  r  >  0,  we  will  denote  by  Xr  the  linear  space  of  all  two-sided  complex 
sequences  a  =  {aTO},  m  —  0, ±1, ±2,  ■  •  •  such  that  for  some  c  >  0, 

l«-l '  S"  •  %/S  <  c  (3.6) 

CT 

for  all  m  >  r.  We  will  denote  by  Yr  the  linear  space  of  all  complex  sequences  0  =  {/?m},  m  = 
0,  ±1,  ±2,  •  •  •  ,  such  that  for  some  c  >  0, 

IA»|(^)’VS<c  (3.7) 

for  all  m  >  r.  It  is  easy  to  see  that  Xr  C  Yr,  and  that  the  condition  (3.6)  is  a  very  restrictive  one, 
since  in  order  to  satisfy  it,  the  elements  of  the  sequence  {om}  must  decay  roughly  as  (r/2)m  •  ml, 
while  the  condition  (3.7)  is  a  very  mild  one  -  it  prohibits  the  elements  of  {/9m}  from  growing 
faster  than  approximately  (2/r)m  •  ml.  By  applying  formulae  (9.3.1),  (9.3.2)  from  [l],  it  is  easy 
to  show  that  in  (3.1),  (3.2),  o  €  Vj*|*  and  0  €  Conversely,  for  any  sequence  a  € 

the  expansion  (3.1)  converges  inside  D,  and  for  any  0  €  Vjjt|ft«  the  expansion  (3.2)  converges 
outside  D.  For  a  natural  n,  we  will  denote  by  T„  a  linear  mapping  Yr  — *  Yr  converting  a  sequence 
a  =  {am},  m  =  0,±1,±2, •••  into  a  sequence  d  =  {dm},  m  =  0,±1,±2, •••  defined  by  the 
formulae 

dm  =  otm  for  |m|  <  n 
dm  =  0  for  |m|>n+l. 

Clearly,  Tn{Yr)  c  A'r,  and  for  obvious  reasons,  we  will  refer  to  Tn  as  truncation. 

For  the  remainder  of  this  section,  Di,  Di,  Dz  will  denote  three  disks  in  R2  such  that  Dz  C  Di 
and  D\  n  Dz  —  0.  The  centers  and  radii  of  these  disks  will  be  denoted  by  cj,  C2,  cz  and  Ri,  R2,  Rz 
respectively.  We  will  denote  the  distance  ||c2  -  ci||  by  P12,  and  the  distance  ||c3  —  ci||  by  piz-  The 
angle  between  the  vector  C2  -  Cj  and  the  x  axis  will  be  denoted  by  #12.  and  the  the  angle  between 
c3  -  ci  and  the  x  axis  will  be  denoted  by  £13.  For  a  point  x  €  R 2,  we  will  denote  ||ar  -  ci||  by  pi, 
(jar  —  C2H  by  P2-  and  |(ar  -  C3IJ  by  pz-  Finally,  the  angles  between  the  vectors  x  -  c\ ,  x  —  C2  and  x  -  c3 
and  the  x  axis  will  be  denoted  by  6\,  62 .  03  respectively. 

Suppose  that  n  is  a  natural  number  or  oc.  We  will  define  a  linear  operator  U"2Cl  :  A'|*.|Rj  — ► 
A*|*|«,  as  follows.  If  A  =  {am},  m  =  0,±1,±2, ••  ■  is  an  element  of  A"|*|/js  then  C^c"ei(A)  =  B  = 
|6m  }i  m  —  0,±1,±2,-  •  ■  is  defined  by  the  formula 

6m  =  e  i)fl3am-jJj{kp\2) 


(3.8) 
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for  all  m  =  0,±1,±2,  •  •  ••  An  operator  Vc"e2  :  “ 1 "  wiU  be  defined  by  the  formula 

bm  =  J2  *-is(tl,-,)am-jMkpn )  (3-9) 

jm-n 

for  any  A  =  { am },  B  =  {bm}  such  that  B  =  V"Cj(A).  Finally,  the  operator  X\k\Rl  —  Vjfc|/ea 
converting  an  arbitrary  sequence  A  =  {am }  €  A')*!*,  into  the  sequence  B  =  {6m  }  €  defined 
by  the  formula 


bm  =  £  e-0'(#i*-')am-^y(fcpi3) 


(3.10) 


will  be  denoted  by  lVe"lC8.  For  any  natural  m,  we  will  define  the  operators  :  A'm«s  — *  A)*^, , 
V™  :  Vj*,*  -  Ym„  W™  :  Xm,  -  Ym,  by  the  formulae 

U™=TmU?3tlTm, 

r»m  _  <T»  .  yn  .  T 

eiej  *  m  ei*2  1  m  1 

ii/nm  _  7.  .  ty"  .  t< 

neie»  —  M  ei«* 

Remark  3.2 

We  will  denote  by  Y  the  set  of  all  two-sided  complex  sequences  {y,},  i  =  ±1,±2,  •••. 
When  n  is  a  natural  number  (as  opposed  to  oo),  the  formulae  (3.8),  (3.9),  (3.10)  define  mappings 
••  Y  -*■  Y  that  can  be  viewed  as  extensions  of  the  mappings  U?3Cl,V£e3,W"et  i.e., 


U”  =17" 

eaCllx|i|Rj  u'* 


V"  —  V"  IV"  —  IV" 

eaci’  eiealA'itm,  ejei»  eic*|X|*|Rl  cits’ 


Whenever  there  is  no  danger  of  confusion,  we  will  make  no  distinction  between  the  mappings 
Kn,ca,  ^c"c,  and  the  mappings  Uen3Cl,  VenjC3,  W?lCt  respectively. 

The  following  three  lemmas  justify  our  referring  to  the  operators  U?3Cl ,  K"eS’  U™\  > 

Ve"” ,  IV,."”  as  translation  operators,  and  will  be  used  to  shift  the  origins  of  H  and  7  -  expansions, 

and  to  convert  77-expansions  into  7-expansions.  They  are  a  direct  consequence  of  the  Graf’s 
addition  theorem  for  Bessel  functions  (see  [l],  (9.1.79)). 

Lemma  3.2.  Suppose  that  ip  :  R2  —*  C1  is  a  radiation  held  analytical  in  R2  \  D2  and  satisfying 
the  radiation  condition  (2.10)  at  oo.  Suppose  further  that  ip  is  represented  by  an  expansion 


ip{x)=  £  PmHm(kp2)eimh 


(3.11) 


valid  in  R2  \  D2,  and  by  an  expansion 


iP(x)=  £  dmHm(kpi)eim9' 


(3.12) 


valid  in  R2\D1.  Then  {8m}  =  U%Cl({/3m}). 

Lemma  3.3.  Suppose  that  ip  :  R2  — •  C1  is  a  radiation  held  analytical  in  D\  and  represented  by  an 


expansion 


*(*)-  £  oimJm{kp\)eimil 


(3.13) 
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valid  in  Di ,  and  by  an  expansion 


+00 

*(*)«  ^  &mJm(kp2)eimt>  (3.14) 

mas— 00 

valid  in  D2.  Then  {am}  =  ^“2({am}). 

Lemma  3.4.  Suppose  that  rp  :  J?2  — +  C1  is  a  radiation  field  analytical  outside  the  disk  Di  and 
satisfying  tie  radiation  condition  (2.10),  and  that  it  is  represented  in  R2  \  D\  by  the  expansion 
(3.12).  Then  inside  the  disk  Da,  tie  function  ip  can  be  represented  in  tiie  form 

+00 

V>(*)  =  E  ~1mJm(kpl3,)e'mta  (3.15) 

m»— oo 


With  hm}  =  W^a&n}). 

Under  the  conditions  of  Lemma  3.2,  we  will  define  a  radiation  field  <p"aej  :  It2  \  D\  — +  C1  by 
the  expression 

E  (3.16) 

m*-oo 

with  the  coefficients  m  =  0,±1,±2,  •  •  •  defined  by  the  formula  {6  m  }  =  U^cxiPm})-  Simi¬ 

larly,  under  the  conditions  of  Lemma  3.3,  we  will  define  a  radiation  field  ip?lC2  :  D2  —*  C1  by  the 
formula 

+00 

«„,<*)  =  £  Sm •WW’**  (3.17) 

m=-oo 

with  {6m}  =  Ve"C2({am}).  Finally,  under  the  conditions  of  Lemma  3.4,  we  will  denote  by  $"lC#  the 
radiation  field  Z>3  — *  C1  defined  by  the  formula 

+00 

&",.,(*)=  E  ^Jm(kps)eimta  (3.18) 


with  the  sequence  (£m}i  m  =  0,±1,±2,-  •  •  defined  by  the  formula  {6m}  =  U7enlCll({#m})- 

Obviously,  (x)  =  ip(x)  for  any  x  €  R2  \  D\,  <pfaCl  (x)  =  <p(x)  for  any  x  €  Da,  and 
tf-,  (x)  =  VK*)  for  any  a:  €  Z?3,  and  we  will  view  the  mappings  0?2Cl,  <A?lCs>  “  approximations 
to  the  mappings  ip,  <t>  and  ip  respectively. 

3.3.  Diagonalization  of  translation  operators. 

Theorem  3.1.  For  any  natural  n,  each  of  the  operators  U"aCl,  VC"C2 ,  H7e" Cg ,  is  a  bounded  scalar 
type  operator  (see  [8],  v.  JIT).  Their  sets  of  eigenvectors  coincide,  and  are  given  by  the  formula 

e,  =  {e'*m},  m  —  0,±1,±2,  •  •  ■  (3.19) 


with  q  €  [0,  2tt] .  The  spectral  representations  of  the  operators  U"sCl,  V'C"C2 ,  H’"lCs  are  respectively 

^e"Cl=  r  K(q)Pqdq, 

Jo 


(3.20) 
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f 2» 

^e"e j  ~  /  PnWjPqdq,  (3.21) 

Jo 

W?,e,=  r^n(q)Pvdq,  (3.22) 

Jo 

where  Pq  is  a  rani  one  operator  projecting  Y  on  its  subspa.ce  spanned  by  the  vector  eq  and  com¬ 
muting  with  each  of  the  operators  I7£ei,  VC"C2,  and  the  functions  An,  pn,  vn  :  [0,27r]  — *  C1 

are  defined  by  the  formulae 


n 


An(?)  =  £  e~,m(*+#**) Jm  (kpi2) , 

m«-n 

(3.23) 

Mn (?)  =  Y,  i2), 

(3.24) 

*n(?)  =  t  c_,m(,+,1*_,r) Hm {kp\z) . 

(3.25) 

Proo/.  An  inspection  of  the  formulae  (3.8),  (3.9),  (3.10)  shows  that  the  mappings  t/£Cl,  VC”CJ, 
are  convolutions  of  the  sequence  {ame},  m  =  0,±1,±2,  •  •  • 

with  the  finite  sequences 

{e-m#”Jm(*p12)},  m  =  0,±1,±2,  •  •  • 

(3.26) 

)jn(kp12)},  m  =  0,4:1, ±2,-  •  • 

(3.27) 

{e-im{llt'v)Hm(kp12)},  m  =  0,±1,±2,  •  •  • 

(3.28) 

respectively.  Convolutions  with  finite  sequences  are  diagonalized  by  the  Fourier  Transformation 
(see,  for  example,  [8,  14]),  which  proves  a).  We  obtain  b),  c),  and  d)  by  applying  the  Fourier 
Transformation  to  the  sequences  (3.26),  (3.27),  (3.28). 


Remark  3.3 

For  the  mappings  I7£Cl,  Ve”ea,  the  above  theorem  can  be  extended  to  the  case  n  =  oo.  As  is 
well  known,  for  any  z  €  Cl,  0  €  [0,27t], 

+00 

J2  =  eitcot*  (3.29) 

m*c— oo 

(see,  for  example,  [l],  9.1.21),  and  for  n  =  oo  the  expressions  (3.20),  (3.21)  assume  the  form 

Aoo(tf)  =  e''*'>2COS<«+#»s\ 

px(q)  =  e,'*',12COs(9+*,2-,r). 

However,  Theorem  3.1  can  not  be  extended  to  the  case  of  the  operator  since  the  latter  is 

unbounded  (see  [8], VIII). 

Remark  3.4 


(3.30) 

(3.31) 
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In  numerical  calculations,  the  operators  C/ ,  with  sufficiently  large  n,  m  will  be 
viewed  as  approximations  to  the  operators  U"aCl ,  V"iea ,  IV” Cj  (see  Remark  3.1).  Clearly,  applying 
either  of  the  operators  I/*®,  Vf“™,  IV”™  numerically  to  an  arbitrary  sequence  is  an  order  nm 
procedure,  which  can  be  prohibitively  expensive  for  large  scale  problems.  However,  it  follows  from 
the  above  theorem  that  the  operators  I/e”™  ,  Ve*J* ,  IV£™  are  convolutions  of  sequences  of  lengths 
n,  m,  and  such  convolutions  can  be  evaluated  in  order  (n  +  m)  log(n  +  m)  operations  by  means  of 
the  Fast  Fourier  Transformation  (see  [4,  12]). 

3.4.  Asymptotic  forms  of  radiation  fields. 

In  this  subsection,  we  introduce  an  alternative  form  of  the  expansions  (3.1),  (3.2)  and  the 
mappings  U^Cl ,  V"es,  W™e»>  providing  a  simple  physical  interpretation  of  the  Lemmas  3.2  -  3.4 
and  the  Theorem  3.1. 

For  the  expansion  (3.2),  we  will  consider  a  function  F^Xo  :  [0,2tt]  — ♦  C1  defined  by  the  formula 

"  ,5®  V>(*-*  +  xo)  ■  >At)  ■  e~ikt  •  •  e«  *  (3.32) 

with  a:=  (cos0,sin0).  Substituting  (2.18),  (3.2),  into  (3.32),  we  obtain 

+oo 

4I0(«)=  £  Aae-V'V^  (3.33) 

m=— oo 

which  provides  an  explicit  expression  of  via  the  coefficients  {/3m},  and  we  will  refer  to  F^,*0 
as  the  asymptotic  representation  of  the  field  0  with  the  origin  at  xo- 

In  order  to  define  an  asymptotic  representation  of  the  field  4>  in  (3.1),  we  will  have  to  introduce 
an  additional  assumption  that 

+oe 

£  |orm|  *  c  <  oc  (3.34) 

m*— oo 

which  guarantees  that  the  expansion  (3.1)  converges  on  the  whole  R2.  By  combining  (2.17)  and 
(3.1),  it  is  easy  to  show  that  for  any  6  €  [0,  2tt]  and  x  =  (cos  6.  sin  6),  there  exist  unique  numbers 
Ug,  Vg  such  that 

lim  0(*  •  x  +  x0)  •  y/(2 fart)  -  (e<*‘-T)’  ■  U,  +  e .  V,  eikt  =  0),  (3.35) 

t—+  00 

and  that  the  numbers  Ug,  Vg  are  defined  by  the  formulae 

+00 

U,=  £  (Om-e-^)-^,  (3.36) 

oo 

-f-co 

Vg  =  £  (am  •  e2?')  e,m# .  (3.37). 

oo 

Now.  we  will  define  the  function  Gc,Xo  :  [0.2ff]  —  C1  by  the  formula 

GW*)  =  CT.  (3.38) 

and  refer  to  it  as  the  asymptotic  representation  of  the  field  0  with  the  origin  at  xo- 

The  following  three  lemmas  establish  that  the  asymptotic  representations  of  the  fields  0,  <t> 
diagonalize  the  translation  operators  17C"CI,  ^C"CJ,  H*"*,.  All  three  are  an  immediate  consequence 
of  the  Lemmas  3.2-3.4  and  the  Theorem  3.1. 
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Lemma  S.5. 

Under  the  conditions  of  Lemma.  3.2,  for  any  8  €  [0,2tt]. 

F \l>,c\{6)  =  F i>,ci(0) '  ^oo (8)  (3.39) 

with  A*,  :  [0,2jt]  —*  C1  defined  by  (3.30). 

Lemma  3.6.  Under  the  conditions  of  Lemma  3.3 ,  for  any  8  6  [0,2?r], 

Gy,c2(0)  =  Gt/\C,(0)  •  Moo(0)  (3.40) 

with  Hoo  ■  [0,2r]  — *  C1  defined  by  (3.31). 

Lemma  3.7.  Under  the  conditions  of  Lemma  3.4,  for  any  e  >  0  there  exists  N  >  0  such  that  for 
any  n  >  N, 

\v(x)  -  4>n(x)\  <  e  (3.41) 

for  any  x  €  D$  where  <j>„  :  R2  — *  C1  is  a  radiation  held  analytical  inside  D$  with  the  asymptotic 
representation  defined  by  the  formula 


Gtt>n,  c»  —  t'n{8)F^,Xl(8),  (3-42) 

and  the  function  vn  :  [0, 27r]  — »  C1  is  defined  by  (3.22). 

The  following  two  lemmas  are  a  direct  consequence  of  the  Theorem  3.1,  Lemmas  3.5,  3.6 
and  the  Remark  3.3.  They  provide  explicit  expressions  for  the  asymptotic  representations  of  the 
fields  generated  by  a  charge  and  a  dipole  located  at  an  arbitrary  point,  and  for  the  values  and  the 
derivatives  of  a  field  given  by  its  asymptotic  representation. 

Lemma  3.8.  Suppose  that  under  the  conditions  of  the  Lemma  3.2,  the  held  ip  is  defined  by  the 
formula  ip  =  <pkx,  i.e.,  it  is  the  held  of  a  unity  charge  located  at  the  point  x.  Then  the  asymptotic 
representation  of  the  held  ip  with  the  center  ci  is  given  by  the  formula 

F<,iei(0)  =  e~ikfiieos{,~fi)-  (3-43) 

If  the  held  ip  is  given  by  the  formula  ip  —  4>kh,  i.e.,  it  is  the  held  of  a  unity  dipole  located  at 
the  point  x  and  oriented  in  the  direction  h,  then  the  asymptotic  representation  of  (3.12)  with  the 
center  at  c\  is  given  by  the  formula 

/V.e,  W  =  •*  cos (8  -  8h)e~ikf ‘  “»(«"*»>  (3.44) 

where  8/,  is  the  angle  between  the  vector  h  and  the  x  axis. 

Lemma  3.9.  Suppose  that  under  the  conditions  of  the  Lemma  3.3,  the  held  (3.13)  has  the  asymp¬ 
totic  representation  G^x ,  :  [0,2"]  —  C1.  Then  for  any  x  €  D\.  h  €  i?2  such  that  ||/i||  =  1. 

4>(x)  =  f<W*)  •  e-^'^-^dO.  (3.45) 

-57T  Jo 

jt4>(x  +  th) |(_0  =  £  fj  (8)  cos (8  -  8h)e-'k<"  ™V~'>)d8.  (3.46) 


3.5.  Numerical  evaluation  of  translation  operators. 

For  the  rest  of  this  paper,  we  will  view  the  asymptotic  representations  (3.32),  (3.38)  (as  opposed 
to  the  expansions  of  the  forms  (3.1),  (3.2))  as  our  principal  tool  for  representing  radiation  fields. 
Lemma  3.8  permits  one  to  calculate  asymptotic  representations  of  fields  of  distributions  of  charges 
and  dipoles  without  evaluating  the  coefficients  of  their  if -expansions,  and  Lemma  3.9  provides  a 
tool  for  calculating  the  fields  and  derivatives  of  the  fields  with  given  asymptotic  representations 
without  having  to  evaluate  the  coefficients  of  J -expansions  of  these  fields. 

For  a  radiation  field  ip  :  R2  — ♦  C1  analytical  outside  Di  and  satisfying  the  radiation  condition 
(2.10),  and  an  integer  n  >  2,  we  will  denote  by  F$  the  finite  sequence  01,02,  •  •  ,On  defined  by 
the  formulae 

a,'  =  (®i)>  (3*47) 

W{  =  2ir  —  —  — ,  t  =  1,2,  •••,«.  (3.48) 

Similarly,  for  a  radiation  field  <f>  analytical  inside  D\  and  possessing  an  asymptotic  representation 
we  will  denote  by  Gjjei  the  sequence  i>i,  62,  •  ■  • ,  bn  defined  by  the  formula 

b»  =  G  (u.',).  (3.49) 

Obviously,  ,  GJ  ei  are  tabulations  at  n  equispaced  nodes  on  the  interval  [0,2ir]  of  the  functions 
/V,e, ,  G^d  respectively,  and  we  will  view  F^Ci ,  GJ  as  finite-dimensional  projections  of  the 
asymptotic  representations  of  the  ip,  <j>  radiation  fields. 

For  a  finite  sequence  GJei  =  {fci ,  f>2.  ■  •  • ,  b„  }  we  will  consider  a  radiation  field  Gj  c  :  R2  -*  C 1 
defined  by  the  formula 

<5;«(i)  =  (3.50) 

Clearly,  (3.50)  is  a  trapezoidal  approximation  to  the  integral  (3.45),  and  we  will  look  upon  (3.50) 
as  an  approximation  to  the  field  <p.  Differentiating  (3.50)  with  respect  to  x,  we  obtain  the  formula 

+  «0|.-o  =  E  5,  cos(wj  -  **)«-“"  “<*/-*■)  (3.51) 

j- 1 

for  any  h  €  R2,  and  we  will  view  (3.51)  as  an  approximation  to  (3.46). 

Finally,  we  will  define  mappings  Pc”",  :  Cn  — +  C"  by  the  formulae 

Fc^ci  '  ' '  3  ^n)  =  (^m(u,l)  '  2i,  ^m(^ 2)  '  22i  ‘  '  3  ^m(Wn)  '  2n),  (3.52) 

Q?£AZ  13*23-* -3^)  =  (Mm(wi)  -2l,Mm(«>2)  •22,--3Mm(«>n)  ‘  *n)5  (3.53) 

5££(*1»«2.‘  •  •  3«n)  =  •  Zi,  (/»(«*)  ■  *2,-  '  *  t^mK.)  •  2»),  (3.54) 

with  the  functions  Am,  nm,  um  defined  by  (3.20),  (3.21),  (3.22).  It  is  easy  to  see  from  the  formulae 
(3.22)-(3.23),  (3.52)-(3.54)  and  the  Theorem  3.1  that  under  the  conditions  of  the  Lemmas  3.2-3.4, 

P>mn  (  jpn  \  __  rn 
C2ClV-*^,C2/  *UiC  j’ 

Q?,ncAGlc,)  =  G?'C2. 

5mn  (  r”  \  _ 

cjca 

with  U  =  ip%Cl,  V  =  <f>dc 2 and  we  w511  look  uP°n  the  operators  P™ ,  S™  as 

discretizations  of  diagonal  forms  of  the  operators  ,  Fe"‘c",  IV”™. 

In  order  to  estimate  the  number  of  nodes  in  the  discretizations  (3.47),  (3.49)  required  to  obtain 
a  given  accuracy  £  in  the  evaluation  of  the  fields  <p,  ip  we  will  need  the  following  well  known  fact 
(see,  for  example  [7]). 
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Lemma  3.10.  For  a ay  integer  m,  n  such  that  n  >  2|m|,  the  n-point  trapezoidal  quadrature  rule 
on  the  interval  [0, 2ir]  integrates  the  function  eimz  exactly. 

Remark  3.5 

By  combining  the  above  lemma  with  Remark  3.1,  it  is  easy  see  that  the  number  n  of  nodes 
in  the  discretization  GJej  of  the  function  G*iCl  :  [0,2ff]  — *  C1  has  to  be  approximately  equal  to 
2|fc|Ri,  and  is  almost  independent  of  the  accuracy  e  with  which  the  field  4>  is  being  calculated. 

Lemmas  3.5  •  3.7  provide  a  tool  for  shifting  the  origins  of  asymptotic  expansions  of  radiation 
fields,  and  for  converting  asymptotic  representations  of  the  form  (3.32)  into  asymptotic  represen¬ 
tations  of  the  form  (3.38)  for  a  cost  proportional  to  n,  where  n  is  the  number  of  nodes  in  the 
discretization  (3.48)  of  the  interval  [0,2tt].  In  the  following  two  sections,  this  apparatus  is  used  to 
construct  an  algorithm  for  rapid  evaluation  of  integral  operators  of  Section  2. 

4.  Rapid  Evaluation  of  Integral  Operators  of  Section  2 

In  this  section,  we  describe  an  algorithm  for  rapid  evaluation  of  the  field  and  the  normal 
derivative  of  the  field  created  on  a  curve  7  by  charge  and  dipole  distributions  on  that  same  curve. 
For  definitiveness,  we  will  be  discussing  the  evaluation  of  the  field  created  by  a  charge  distribution. 
The  algorithms  evaluating  the  normal  derivative  of  the  field  created  by  a  charge  distribution,  and 
the  field  and  the  normal  derivative  of  the  field  created  by  a  dipole  distribution  are  quite  similar. 

4.1.  Notation 

We  will  consider  the  situation  depicted  in  Figure  2.  The  curve  7  is  discretized  into  equispaced 
nodes  £1,  £2,  •  •  •,  x„,  and  we  will  denote  the  spacing  between  the  adjacent  nodes  by  h.  Suppose 
that  for  each  i  —  1,2,  •  •  • ,  n,  a  charge  a,-  of  strength  cr,  is  located  at  the  point  x,-.  In  this  section, 
we  describe  an  algorithm  for  rapid  calculation  of  approximations  <7,  =  1,2,  •  •  • ,  n  to  the  sums 

n 

G*{xi)  =  (41) 

i-i 

i+i 

for  t  =  1,2.  •  •  •  ,n.  Clearly,  this  is  an  order  n2  process  (evaluating  n  fields  at  n  points).  However, 
if  we  are  interested  in  evaluating  (4.1)  with  a  finite  accuracy  (which  is  always  the  case  in  practical 
calculations),  Theorem  3.1  and  Lemmas  3.7,  3.8,  3.9  can  be  used  to  speed  up  the  process. 

For  an  integer  m  >  2,  we  will  define  the  points  <i,*2>- •  •  ,*m+i  on  the  interval  [0,L]  by  the 
formula  U  —  (t-  l)L/m,  subdividing  the  interval  [0,  L]  into  m  segments  of  equal  length,  and  denote 
the  center  of  the  i-th  segment  by  z,-,  so  that  z,  =  U  +  L/(2m).  For  each  natural  j  =  1,2,  •  •  •  ,m, 
we  will  denote  by  Aj  the  set  of  all  charges  a,-  such  that  x,-  €  7 ([</,  </+i]),  and  by  Dj  the  circle  of 
radius  r  =  L/(2m)  with  the  center  at  7 (zy).  We  will  denote  by  Wj  the  union  of  all  Aj  such  that 
|| Zj  -  z,j|  >  3r,  and  by  Wj  the  union  of  all  Aj  such  that  ||zy  -  z,j|  <  3r.  Obviously,  Aj  C  Dj  for  any 
j  =  1,2,  •  •  •  ,m.  Also,  it  follows  from  the  triangle  inequality  that 

min  ||x  -  y\\  >  r 

z£A, 

y&Aj 

for  any  i,j  such  that  Aj  C  Wj.  Finally,  we  will  denote  by  Oj  the  field  of  all  charges  a,-  such  that 
x,  C  Aj  and  observe  that  if  xp  €  Aj  then 

G<r(Xp)  —  )  ]  ^>{xp)  +  y*!  °V ^Xjfep)’  (4*2) 

A,zWj  *,€»•> 
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4.2.  Detailed  description  of  an  order  n3/2  algorithm. 

In  this  subsection,  M,  N  will  denote  “sufficiently  large*  integer  numbers.  The  actual  choice  of 
the  numbers  M,  N  is  discussed  in  the  following  subsection. 

We  will  evaluate  the  fields  (4.1)  in  five  steps. 

Step  1. 

Using  Lemma  3.8,  obtain  discretized  asymptotic  representations  F£  of  the  fields  <f>j  for 
all  j  =  1,2,-  •  -,m. 

Step  2 

For  every  pair  of  natural  numbers  i,j  €  [l,m]  such  that  A,-  C  Wj,  calculate  the  representation 


pA’  _  g>M,N  i  pN  \ 


of  the  field  faj  =  j  and  view  it  as  a  finite-  dimensional  approximation  to  the  asymptotic 

representation  of  the  field  fa  on  Dj. 

Step  3 

For  each  natural  j  €  [l,m],  calculate  the  sum 


AiCWj 


and  view  the  field  ipj  =  Yli'l’ij  88  811  approximation  to  the  field  and  as  a 

finite-dimensional  approximation  to  the  asymptotic  representation  of  rpj  on  Dj. 

Step  4 

For  each  natural  j  €  [l,m],  evaluate 


for  all  t  such  that  x ,■  G  7([^,*>+i])  and  look  upon  (4.5)  as  an  approximation  to  rpj{xi).  t 
Step  5. 

For  each  j  =  1,2 ,  ■  •  • ,  m,  evaluate  the  sum 


M*i)+  H  °P<t>kz,{xi) 


for  all  i  such  that  x,-  €  'y([f/*f/+i])*  and  view  (4.6)  as  an  approximation  to  G<r(x,). 

4.3.  Choice  of  parameters  and  CPU  time  estimate 

In  the  estimates  below,  0,  b,  c,  d,  e  are  coefficients  determined  by  the  computer  system,  language, 
particular  implementation  of  the  algorithm,  etc. 

Step  1 

Obviously,  this  step  will  require  order  n  ■  A’  operations  (tabulating  F£.  at  N  nodes  on  the 
interval  [0,2?r]  for  each  of  the  nodes  xi,X2,-  -,x„  ).  According  to  the  Remark  3.5,  N  ~  |fc|  •  L/m. 
and  the  CPU  time  estimate  for  this  step  becomes  a  ■  n  •  |fc|  •  L/m. 

Step  2 

For  each  of  the  pairs  i,j  such  that  Aj  C  H’,,  evaluating  (4.3)  will  require  order  N  operations 
(see  (3.54)  ),  and  the  total  number  of  such  pairs  is  less  than  m2,  which  results  in  the  CPU  time 
estimate  of  6  •  m2  •  n  ~  b  ■  m2  •  \k\  •  L/m  =  b  •  m  •  (&|  •  L  for  this  step. 

Step  3 


■SAW 


SUM 


Obviously,  evaluating  the  sums  (4.4)  for  all  j  —  1,2,  ••  •  ,m  is  an  order  cm-N  ~  cm\k\L/m  = 
c-\k\-  L  procedure. 

Step  4 

Evaluating  (4.5)  for  each  »  =  1,2,  •  •  • ,  n  is  an  order  N  procedure,  resulting  in  the  total  CPU 
time  estimate  for  this  step  of  d  ■  n  ■  N  ~  d  ■  n  ■  |fc|  •  L/m. 

Step  5 

Evaluating  the  sum  (4.6)  for  each  t  =  1,2,  •  •  •  ,n  is  an  order  n/m  procedure,  with  the  resulting 
CPU  time  estimate  of  t  ■  n2/m  for  this  step. 

Summing  up  the  time  estimates  for  the  steps  1-5,  we  obtain  the  following  time  estimate  for 
the  whole  process: 


T  =  A  ■  n  ■  }Ar|  •  L/m  +  b  ■  m  ■  |Jt|  •  L  +  c  ■  |Jt|  •  L  + 


e  •  nr 
m 


(4.7) 


with  A  =  fl  +  i,  and  we  would  like  to  choose  m  in  such  a  manner  that  (4.7)  would  be  minimized. 
Differentiating  (4.7)  with  respect  to  m,  and  setting  the  resulting  derivative  to  zero,  we  obtain 


A  ■  n  ■  Ifcl  •  L  +  e  •  n2 
mmin  “  V  b  ■  |ik|  •  L 

and  the  corresponding  minimum  of  (4.7)  is  equal  to 

TnUn  =  2 s/a  •  n  •  |Jfe|  •  I  +  e  •  »a  •  y/b-\k\-L+c  ■  k  ■  L. 


(4.8) 


(4.9) 


If  the  calculations  are  performed  with  a  fixed  number  of  nodes  per  wavelength  (which  is  often 
a  reasonable  assumption),  n  is  proportional  to  |&|£,  and  (4.9)  assumes  the  form 


7n>m  ~  (k  ■  L)$ 


(4.10) 


or 

rroi„  ~  n§  (4.11) 

which  for  large  n  is  considerably  smaller  than  n2. 

4.4.  Further  reduction  of  the  CPU  time  estimate  of  the  process 

The  approach  of  the  above  subsection  can  be  used  recursively  by  subdividing  each  of  the  sets 
Ai  into  subsets  {B,/},  j  =  1, 2,  •  •  • ,  fh  with  appropriately  chosen  fh  and  representing  the  fields  <j>i  as 
sums  fa  —  4>ij  where  4>ij  is  the  field  created  by  all  charges  ap  such  that  ap  €  B,y.  A  calculation 
similar  to  the  one  in  the  preceeding  section  shows  that  such  an  algorithm  will  have  an  asymptotic 
CPU  time  estimate  of  n4/3. 

In  [18],  such  subdivision  process  is  used  recursively  to  obtain  an  order  n  algorithm  for  numerical 
evaluation  of  integral  operators  of  classical  potential  theory  (Laplace's  equation).  By  reproducing 
the  construction  of  Section  VII  of  [18]  almost  literally,  one  can  obtain  an  order  n  log(n)  algorithm  for 
evaluating  (4.1).  However,  our  estimates  indicate  that  for  problems  of  practicable  size  (n  <  20,000), 
the  improvement  in  actual  computation  times  obtained  by  replacing  an  order  n 4/3  algorithm  with 
an  order  n  log  n  algorithm  would  not  be  very  significant. 

5.  Implementation  and  Numerical  Results 
5.1.  Implementation. 

In  implementing  the  algorithm  of  the  present  paper,  we  closely  followed  [16],  as  far  as  the  for¬ 
mulation  of  the  system  of  integral  equations  (2.13),  (2.14)  and  their  discretization  are  concerned. 
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Equations  (2.13),  (2.14)  were  discretized  by  means  of  the  Nystrom  algorithm,  which  calls  for  re¬ 
placing  the  integrals  with  finite  quadrature  formulae  and  the  integral  equations  with  systems  of 
linear  algebraic  equations  respectively,  (see  [3,  16]).  Since  the  kernels  of  the  equations  (2.13),  (2.14) 
are  logarithmically  singular,  quadrature  formulae  suitable  for  functions  with  logarithmic  singular¬ 
ities  have  to  be  employed,  and  we  used  a  fourth  order  convergent  formula  of  the  type  discussed 
in  [id])-  Discretization  of  the  equations  (2.13),  (2.14)  by  means  of  the  Nystrom  algorithm  leads 
to  large  scale,  asymptotically  well-conditioned  systems  of  linear  algebraic  equations  (see  [3]),  to 
which  we  apply  a  version  of  the  Generalized  Conjugate  Residual  algorithm  (see  [11]).  Solution  of 
a  linear  system  by  means  of  a  Generalized  Conjugate  Residual  type  algorithm  involves  application 
of  the  matrix  of  the  system  to  a  sequence  of  recursively  generated  vectors,  which  is  accomplished 
by  means  of  the  algorithm  of  the  preceeding  section  (we  implemented  an  order  n4/3  version  of  the 
process).  As  is  well  known,  in  order  to  converge  to  a  relative  accuracy  e,  a  Generalized  Conjugate 
Residual  type  algorithm  requires  order  —K  •  log(e)  iterations,  where  K  is  the  condition  number  of 
the  matrix  of  the  system  being  solved.  This  results  in  the  total  CPU  time  estimate  of  the  order 
n!  •  log(e)  for  the  solution  of  the  equations  (2.13),  (2.14),  since  the  application  of  the  Nystrom 
algorithm  to  these  equations  leads  to  well-conditioned  linear  systems  (see  [3]).  After  the  system 
of  equations  (2.13),  (2.14)  is  solved,  evaluating  the  scattered  field  at  any  point  x  T  involves  two 
numerical  integrations  over  T  (evaluating  the  field  of  the  charge  distribution  and  that  of  the  dipole 
distribution),  which  is  an  order  n  procedure.  This  brings  us  to  the  estimate 

Tsolve  —  ‘  +  b  -  n  ■  m 

for  obtaining  the  solution  of  the  problem  (2.6)  -  (2.10),  where  m  is  the  number  of  points  in  R 2 
where  the  solution  has  to  be  calculated. 

S.2.  Numerical  Results. 

We  have  applied  the  algorithm  of  the  present  paper  to  several  acoustic  scattering  problems  in 
two  dimensions.  Three  examples  are  presented  below. 

a.  Scattering  from  a  circular  inclusion.  In  this  case,  the  problem  possesses  an  analytical 
solution,  providing  a  convenient  way  to  verify  the  accuracy  of  the  algorithm.  In  our  example,  the 
radius  of  the  scatterer  was  100  ft.,  and  its  center  was  at  the  origin.  The  densities  inside  and  outside 
the  scatterer  were  1.2  and  1.0  respectively,  and  the  speeds  of  sound  were  20  and  18  respectively.  The 
ambient  field  was  generated  by  a  cylindrically  symmetric  source  located  at  the  point  (100,100),  and 
the  scattered  field  obtained  by  the  algorithm  of  the  present  paper  was  compared  to  the  analytical 
result  at  40  equispaced  points  located  on  the  circle  of  radius  110  with  the  center  at  the  origin  (see 
Figure  3).  The  calculation  was  performed  with  the  angular  frequency  of  the  source  varying  from  10 
to  320,  and  in  all  cases  the  boundary  of  the  scatterer  was  discretized  at  10  nodes  per  wavelength, 
which  usually  results  in  three  to  four  digit  accuracy.  In  Table  1,  u  is  the  angular  frequency,  n  is  the 
number  of  nodes  on  the  boundary  of  the  scatterer  (at  10  nodes  per  wavelength),  e  is  the  resulting 
mean  square  error  at  40  receivers,  and  m  is  the  number  of  iterations  the  Generalized  Conjugate 
Residual  algorithm  took  to  converge  to  4-digit  accuracy.  Following  observations  can  be  made  from 
Table  1. 

1..  The  accuracy  of  the  solution  obtained  with  a  given  number  of  nodes  per  wavelength  is  virtually 
independent  of  the  frequency. 

2..  The  number  of  iterations  required  by  the  GCRA  to  obtain  a  given  accuracy  is  almost  indepen¬ 
dent  of  the  frequency. 

3..  The  actual  CPU  times  required  by  the  algorithm  grow  somewhat  eratically  as  a  function  of 
the  number  of  nodes  on  the  boundary,  but  seem  to  be  in  agreement  with  the  theoretical  CPU 
time  estimate  of  the  algorithm. 
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b.  Scattering  from  a  lense-shaped  object.  In  the  situation  depicted  in  Figure  4,  the  speed  of 
sound  in  the  containing  space  is  30  ft/sec,  and  the  speed  of  sound  inside  the  scatterer  is  20  ft/sec. 
Both  upper  and  lower  surfaces  of  the  scatterer  are  circular  with  the  radius  of  curvature  1  ft,  so 
that  the  focal  distance  of  the  resulting  lense  is  2  feet.  Three  cylindrically  symmetric  sources  are 
located  at  the  points  I,  II,  and  III  respectively  have  the  frequency  F  =  7162  Hz,  resulting  in  the 
wavelengths  0.0041888  outside  the  scatterer,  and  0.0027925 ft  inside  the  scatterer.  The  amplitude 
of  the  field  generated  by  this  configuration  on  the  screen  is  depicted  in  Figure  5.  Clearly,  the  laws  of 
geometrical  optics  should  be  applicable  to  this  situation  with  a  reasonable  degree  of  accuracy,  since 
the  lense  is  about  240  exterior  (or  360  interior)  wavelengths  in  diameter,  and  a  careful  examination 
of  the  Figure  5  shows  this  to  be  the  case.  In  order  to  produce  the  Figure  5,  the  surface  of  the  lense 
in  Figure  4  was  discretized  into  6134  nodes  (with  10  nodes  per  wavelength),  resulting  in  a  system 
of  linear  algebraic  equations  of  dimension  12268.  The  solution  of  this  problem  took  3358  seconds 
of  CPU  time  on  a  VAX-8600. 

c.  Radiation  from  a  worm-shaped  object.  In  this  numerical  experiment,  a  source  of  angular 
frequency'  w  =  640  was  located  inside  an  inclusion  of  a  complicated  shape  in  a  fluid  two-dimensional 
space  (see  Figure  6  ).  The  densities  inside  and  outside  the  inclusion  were  1.25  and  1.0  respectively, 
and  the  speeds  of  sound  were  20  and  25  respectively.  The  boundary  of  the  inclusion  was  discretized 
at  10  nodes  per  wavelength,  resulting  in  the  total  number  of  5078  nodes  on  the  boundary  of  the 
inclusion.  The  calculation  took  3867  seconds  of  CPU  time,  and  the  amplitude  of  the  field  as 
measured  by  receivers  in  Figure  is  displayed  in  Figure  7. 

6.  Generalisations  and  Conclusions 

6.1.  Other  boundary  conditions. 

In  the  present  paper,  we  have  been  solving  the  problem  (2.6)  -  (2.10),  corresponding  to  the 
scattering  of  sound  from  a  fluid  inclusion  in  a  fluid  space.  Obvously,  the  two  classical  acoustical 
scattering  problems  (Dirichlet  and  Neumann  problems  for  the  Helmholtz  equation,  corresponding  to 
scattering  from  a  cavity  and  from  a  rigid  inclusion  respectively)  can  be  solved  in  a  similar  manner. 
Generally,  whenever  a  scattering  problem  is  reduced  to  a  set  of  second  kind  integral  equations  on 
the  boundary  of  the  scatterer  by  means  of  some  form  of  Stoke’s  theorem  ,  the  algorithm  of  the 
Section  IV  provides  a  tool  for  solving  these  equations  in  order  n4/3  operations. 

6.2.  Multiple  scatterers. 

The  case  of  multiple  scatterers  does  not  differ  substantially  from  the  case  of  a  single  scatterer. 
The  scattered  field  inside  each  scatterer  is  represented  in  a  manner  precisely  similar  to  the  one  used 
to  solve  a  single  scatterer  problem.  The  scattered  field  outside  the  scatterers  is  represented  by  the 
sum  of  the  fields  of  single  and  double  layer  potentials  on  the  surfaces  of  all  scatterers.  The  whole 
procedure  is  completely  strightforward. 

6.3.  Three-dimensional  version  of  the  theory. 

Three-dimensional  equivalents  of  the  expansions  (3.1),  (3.2)  are  well  known  (see,  for  example, 
[15]),  and  those  of  the  mappings  U?iey,  Ve"C3 ,  H7cnlCg  are  fairly  easy  to  define.  However,  in  two 
dimensons,  the  translation  operators  UcnaC|,  W’cn)Cg  are  convolutions,  and  the  diagonal  form 

of  the  latter  is  well  known  from  the  standard  theory  of  the  Fourier  integral  (see  Theorem  3.1). 
In  three  dimensions,  no  such  ready-made  analytical  apparatus  is  available,  and  the  proof  of  an 
equivalent  of  the  Theorem  3.1  is  considerably  more  involved.  Fortunately,  diagonal  forms  of  the 
three-dimensional  analogues  of  the  operators  vy]es,  H’”eg  can  be  obtained  by  a  generalization 
of  the  technique  used  in  [20]  to  derive  Graf’s  addition  theorem,  permitting  fast  solvers  for  scattering 
problems  in  three  dimensions  to  be  constructed.  Details  of  this  generalization  will  be  reported  at 
a  later  date. 
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6.4.  Conclusions 

One  of  principal  difficulties  arizing  in  the  solution  of  large-scale  scattering  problems  by  means 
of  integral  equations  is  the  fact  that  the  Green’s  function  for  the  Helmholtz  equation  decays  slowly. 
As  a  result,  the  kernels  of  the  obtained  integral  equations  are  not  sparse,  and  their  discretization 
leads  to  dense  large  scale  systems  of  linear  algebraic  equations.  Solution  of  such  systems  is  an 
extremely  expensive  process  (see  (9,  11]),  which  limits  the  usefulness  of  the  whole  approach. 

In  Section  4  of  the  present  paper,  we  construct  an  algorithm  for  rapid  application  of  the 
matrices  resulting  from  discretization  of  integral  equations  of  scattering  theory  to  arbitrary  vectors. 
The  asymptotic  CPU  time  estimate  of  the  algorithm  is  n4/3,  and  can  be  reduced  to  n  log(n)  where 
n  is  the  number  of  nodes  in  the  discretization  of  the  boundary  of  the  scatterer.  By  combining  the 
approach  of  the  Section  4  with  a  Generalized  Conjugate  Residual  type  process,  we  obtain  an  order 
ns  log(e)  algorithm  for  the  solution  of  the  integral  equations  (2.13),  (2.14)  of  acoustic  scattering 
theory.  This  results  in  acceptable  computation  times  even  for  large-scale  scattering  problems,  (see 
preceeding  section),  as  long  as  the  solution  has  to  be  evaluated  at  a  limited  number  of  points. 
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A  Fluid  Scatterer  Imbedded  in  a  Two-dimensional  Fluid  Space 


Figure  5 


